function p_ri = ri_fe(iter,beta,sigma,y,X,stations,tail)

    station1 = stations(:,1);
    station2 = stations(:,2);
    station3 = stations(:,3);
    
    sim_data = csvread('stations_randsim.csv',1,0);

    stationsim = sim_data(:,1);
    treatsim = sim_data(:,2:10001);

[N,k] = size(X);

treat1 = zeros(N, iter);
treat2 = zeros(N, iter);
treat3 = zeros(N, iter);

treatany = zeros(N,iter);
b = zeros(k+1,iter);
se=zeros(k+1,iter);
test = zeros(k+1,iter);

%Getting ac-wise treatment status
for j = 1:iter
    for i = 1:N
        for s = 1:60 
            if stationsim(s)== station1(i)
                treat1(i,j)=treatsim(s,j);
            end;
            if stationsim(s)==station2(i)
                treat2(i,j)=treatsim(s,j);
            end;
            if stationsim(s)==station3(i)
                treat3(i,j)=treatsim(s,j);         
            end;
            if treat1(i,j)+treat2(i,j)+treat3(i,j)>0
                treatany(i,j)=1;
            end;
        end;
    end;
end;

for a = 1:iter
    [b(:,a),se(:,a)] = fereg(y,treatany(:,a),X,stations,tail);
end; 

%p-value
if tail == 1
    for l = 1:iter
       
            if b(1,l)/se(1,1) < beta(1)/sigma(1)
                test(1,l)=1;
            
            end;
    end;
    p_ri = mean(test,2);
end;

if tail == 2
     for l = 1:iter
        
            if abs(b(1,l)/se(1,l)) > abs(beta(1)/sigma(1))
                test(1,l)=1;
             
            end;
     end;
     p_ri = mean(test,2);
end;